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Abstract — NIFTY, "Numerical Information Field Theory", is 
a software package designed to enable the development of signal 
inference algorithms that operate regardless of the underlying 
spatial grid and its resolution. Its object-oriented framework is 
written in Python, although it accesses libraries written in 
Cython, C++, and C for efficiency. NIFTy offers a toolkit 
that abstracts discretized representations of continuous spaces, 
fields in these spaces, and operators acting on fields into classes. 
Thereby, the correct normalization of operations on fields is 
taken care of automatically without concerning the user. This 
allows for an abstract formulation and programming of inference 
algorithms, including those derived within information field 
theory. Thus, NIFTY permits its user to rapidly prototype 
algorithms in ID, and then apply the developed code in higher- 
dimensional settings of real world problems. The set of spaces 
on which NIFTy operates comprises point sets, n-dimensional 
regular grids, spherical spaces, their harmonic counterparts, 
and product spaces constructed as combinations of those. The 
functionality and diversity of the package is demonstrated by 
a Wiener filter code example that successfully runs without 
modification regardless of the space on which the inference 
problem is defined. 

Index Terms — Inference algorithms, Image reconstruction, In- 
formation theory, Scientific computing, Data analysis, Filtering 
algorithms 

*Email: mselig@mpa-garching . mpg . de 
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I. Introduction 

In many signal inference problems, one tries to reconstruct 
a continuous signal field from a finite set of experimental 
data. The finiteness of data sets is due to their incompleteness, 
resolution, and the sheer duration of the experiment. A further 
complication is the inevitability of experimental noise, which 
can arise from various origins. Numerous methodological 
approaches to such inference problems are known in modern 
information theory founded by Ref. (IH3). 

Signal inference methods are commonly formulated in an 
abstract, mathematical way to be applicable in various scenar- 
ios; i.e., the method itself is independent, or at least partially 
independent, of resolution, geometry, physical size, or even 
dimensionality of the inference problem. It then is up to the 
user to apply the appropriate method correctly to the problem 
at hand. 

The authors are with the Max Planck Institute for Astrophysics (Karl- 
Schwarzschild-StraBe 1, 85748 Garching, Germany). 

Excerpts of this paper are part of the NIFTy source code and documenta- 
tion. 



In practice, signal inference problems are solved numeri- 
cally, rather than analytically. Numerical algorithms should try 
to preserve as much of the universality of the underlying infer- 
ence method as possible, given the limitations of a computer 
environment, so that the code is reuseable. For example, an 
inference algorithm developed in astrophysics that reconstructs 
the photon flux on the sky from high energy photon counts 
might also serve the purpose of reconstructing two- or three- 
dimensional medical images obtained from tomographical 
X-rays. The desire for multi-purpose, problem-independent 
inference algorithms is one motivation for the NIFTY package 
presented here. Another is to facilitate the implementation of 
problem specific algorithms by providing many of the essential 
operations in a convenient way. 

NIFTy stands for "Numerical Information Field Theory". 
It is a software package written in Pytho^ however, it also 
incorporates Cytho^J|4, 5 ], C++, and C libraries for efficient 
computing. 

The purpose of the NIFTY library is to provide a toolkit that 
enables users to implement their algorithms as abstractly as 
they are formulated mathematically. NIFTy's field of applica- 
tion is kept broad and not bound to one specific methodology. 
The implementation of maximum entropy (6l [7), likelihood- 
free, maximum likelihood, or full Bayesian inference methods 
EJE] E) are f eas ible, as well as the implementation of posterior 
sampling procedures based on Markov chain Monte Carlo 
procedures |[TOl [Till . 

Although NIFTY is versatile, the original intention was the 
implementation of inference algorithms that are formulated 
methodically in the language of information field theor)^] 
(IFT). The idea of IFT is to apply information theory to 
the problem of signal field inference, where "field" is the 
physicist's term for a continuous function over a continuous 
space. The recovery of a field that has an infinite number 
of degrees of freedom from finite data can be achieved by 
exploiting the spatial continuity of fields and their internal 
correlation structures. The framework of IFT are detailed in 
Ref. fl2ll where the focus lies on a field theoretical approach to 
inference problems based on Feynman diagrams. An alterna- 
tive approach using entropic matching based on the formalism 
of the Gibbs free energy can be found in Ref. [ 1 3 1 . IFT 



^YTHON homepage http: / /www. python . org/ 

2 CYTHON homepage http://cython.org/ 

3 IFT homepage http: / /www. mpa-garching.mpg.de/ift/ 
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based methods have been developed to reconstruct signal fields 
without a priori knowledge of signal and noise correlation 
structures O [T5l . Furthermore, IFT has been applied to a 
number of problems in astrophysics, namely to recover the 
large scale structure in the cosmic matter distribution using 
galaxy counts fT6U20lL and to reconstruct the Faraday rotation 
of the Milky Way l2T1l . A more abstract application has been 
shown to improve stochastic estimates such as the calculation 
of matrix diagonals by sample averages [ 22 1. 

One natural requirement of signal inference algorithms is 
their independence of the choice of a particular grid and 
a specific resolution, so that the code is easily transferable 
to problems that are similar in terms of the necessary in- 
ference methodology but might differ in terms of geometry 
or dimensionality. In response to this requirement, NIFTY 
comprises several commonly used pixelization schemes and 
their corresponding harmonic bases in an object-oriented 
framework. Furthermore, NIFTY preserves the continuous 
limit by taking care of the correct normalization of operations 
like scalar products, matrix- vector multiplications, and grid 
transformations; i.e., all operations involving position integrals 
over continuous domains. 

The remainder of this paper is structured as follows. In 
Sec. [II] an introduction to signal inference is given, with the 
focus on the representation of continuous information fields 
in the discrete computer environment. Sec. [Ill] provides an 
overview of the class hierarchy and features of the NIFTY 
package. The implementation of a Wiener filter algorithm 



demonstrates the basic functionality of NIFTY in Sec. [TV 
We conclude in Sec. [V] 



II. Concepts of Signal Inference 
A. Fundamental Problem 

Many signal inference problems can be reduced to a single 
model equation, 



d=/(a,...), 



(1) 



where the data set d is the outcome of some function / being 
applied to a set of unknowns]^] Some of the unknowns are 
of interest and form the signal s, whereas the remaining are 
considered as nuisance parameters. The goal of any inference 
algorithm is to obtain an approximation for the signal that is 
"best" supported by the data. Which criteria define this "best" 
is answered differently by different inference methodologies. 

There is in general no chance of a direct inversion of 
Eq. ([T]). Any realistic measurement involves random processes 
summarized as noise and, even for deterministic or noiseless 
measurement processes, the number of degrees of freedom of 
a signal field typically outnumbers those of a finite data set 
measured from it, because the signal field of interest might 
represent a continuous field; e.g., some physical flux or density 
distribution. 

4 An alternative notation commonly found in the literature is y = f[x\. We 
do not use this notation in order to avoid confusion with coordinate variables, 
which in physics are commonly denoted by x and y. 



In order to clarify the concept of measuring a continuous 
signal field, let us consider a linear measurement by some 
response R with additive and signal independent noise n, 



d = Rs + n. 
which reads for the individual data points, 



dx Ri(x)s(x) + rii. 



(2) 



(3) 



Here we introduced the discrete index i G {1, . . . , N} C N 
and the continuous position x G ft of some abstract position 
space ft. For example, in the context of image reconstruction, 
i could label the TV image pixels and x would describe real 
space positions. 

The model given by Eq. ^ already poses a full inference 
problem since it involves an additive random process and a 
non-invertible signal response. As a consequence, there are 
many possible field configurations in the signal phase space 
that could explain a given data set. The approach used to single 
out the "best" estimate of the signal field from the data at 
hand is up to the choice of inference methodology. However, 
the implementation of any derived inference algorithm needs 
a proper discretization scheme for the fields defined on ft. 
Furthermore, it is worthwhile to keep the implementation 
flexible with respect to the characteristics of ft. 

B. Discretized Continuum 

The representation of fields that are mathematically defined 
on a continuous space in a finite computer environment is 
a common necessity. The goal hereby is to preserve the 
continuum limit in the calculus in order to ensure a resolution 
independent discretization. 

Any partition of the continuous position space ft (with 
volume V) into a set of Q disjoint, proper subsets ft q (with 
volumes V q ) defines a pixelization, 



ft 



V 



with q G {1,..., Q} C N, 



dx 



E 



dx 



En 



(4) 



(5) 



Here the number Q characterizes the resolution of the pix- 
elization, and the continuum limit is described by Q — » 00 
and V q — » for all q G {1, . . . , Q} simultaneously. More- 
over, Eq. ([5} defines a discretization of continuous integrals, 
f n dx^J2 q V q . 

Any valid discretization scheme for a field s can be de- 
scribed by a mapping, 



s(x G ft q ) H> S q 



dx W q (x)s(x), 



(6) 



if the weighting function w q (x) is chosen appropriately. In 
order for the discretized version of the field to converge 
to the actual field in the continuum limit, the weighting 
functions need to be normalized in each subset; i.e., Mq : 
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TABLE I 

Overview of derivatives of the NIFTy space class, the corresponding grids, and conjugate space classes. 



NIFTY subclass 


corresponding grid 


conjugate space class 


point_space 


unstructured list of points 


(none) 


rg_space 


n-dimensional regular Euclidean grid over T n 


rg_space 


lm_space 


spherical harmonics 


gl_space or hp_space 


gl_space 


Gauss -Legendre grid on the S 2 sphere 


lm_space 


hp_space 


HEALPix grid on the *S 2 sphere 


lm_space 


nested_space 


(arbitrary product of grids) 


(partial conjugation) 



j Q dx w q (x) = 1. Choosing such a weighting function that 
is constant with respect to x yields 

J n dx s(x) 



k n dx 



(7) 



which corresponds to a discretization of the field by spa- 
tial averaging. Another common and equally valid choice 
is w q (x) = S(x — x q ), which distinguishes some position 
x q Efl q , and evaluates the continuous field at this position, 



/ dx S(x 



X q )s(x) = S(x q ). 



(8) 



In practice, one often makes use of the spatially averaged 



pixel position, 



(x) n ; cf. Eq. ([7]). If the resolution 



is high enough to resolve all features of the signal field s, 
both of these discretization schemes approximate each other, 
(s(x)) n w s((x) n ), since they approximate the continuum 
limit by construction]^] 

All operations involving position integrals can be normal- 
ized in accordance with Eqs. §5§ and ([7]). For example, the 
scalar product between two fields s and u is defined as 



f Q 

S^U = dx S*(x)u(x) ~ 2_) Vq S q U qi 



(9) 



where f denotes adjunction and * complex conjugation. Since 
the approximation in Eq. ^ becomes an equality in the 
continuum limit, the scalar product is independent of the 
pixelization scheme and resolution, if the latter is sufficiently 
high. 

The above line of argumentation analogously applies to the 
discretization of operators. For a linear operator A acting on 
some field s as As = J Q dyA(x, y)s(y), a matrix representa- 
tion discretized in analogy to Eq. ^ is given by 

a, n ^^ a ffn v n a dxd y A (x,y) 
A{x e lip, y e U q ) ^ A pq = 



IIn p n q dxd V 



(10) 



Consequential subtleties regarding operators are addressed in 
App.0 

The proper discretization of spaces, fields, and operators, 
as well as the normalization of position integrals, is essential 
for the conservation of the continuum limit. Their consistent 
implementation in NIFTY allows a pixelization independent 
coding of algorithms. 

5 The approximation of {s(x)) n « s(x q £ Q q ) marks a resolution 
threshold beyond which further refinement of the discretization reveals no 
new features; i.e., no new information content of the field s. 



III. Class and Feature Overview 

The NIFTY library features three main classes: spaces 
that represent certain grids, fields that are defined on spaces, 
and operators that apply to fields. In the following, we will 
introduce the concept of these classes and comment on further 
NIFTY features such as operator probing. 



A. Spaces 

The space class is an abstract class from which all other 
specific space subclasses are derived. Each subclass represents 
a grid type and replaces some of the inherited methods with 
its own methods that are unique to the respective grid. This 
framework ensures an abstract handling of spaces independent 
of the underlying geometrical grid and the grid's resolution. 

An instance of a space subclass represents a geometrical 
space approximated by a specific grid in the computer environ- 
ment. Therefore, each subclass needs to capture all structural 
and dimensional specifics of the grid and all computationally 
relevant quantities such as the data type of associated field 
values. These parameters are stored as properties of an instance 
of the class at its initialization, and they do not need to be 
accessed explicitly by the user thereafter. This prevents the 
writing of grid or resolution dependent code. 

Spatial symmetries of a system can be exploited by cor- 
responding coordinate transformations. Often, transformations 
from one basis to its harmonic counterpart can greatly reduce 
the computational complexity of algorithms. The harmonic 
basis is defined by the eigenbasis of the Laplace operator; 
e.g., for a flat position space it is the Fourier basis]^] This 
conjugation of bases is implemented in NIFTY by distin- 
guishing conjugate space classes, which can be obtained 
by the instance method get_co domain (and checked for 
by check_codomain). Moreover, transformations between 
conjugate spaces are performed automatically if required. 

Thus far, NIFTY has six classes that are derived from the 
abstract space class. These subclasses are described here, and 
an overview can be found in Tab. U 

• The point_space class merely embodies a geometri- 
cally unstructured list of points. This simplest possible 
kind of grid has only one parameter, the total number of 
points. This space is thought to be used as a default data 
space and neither has a conjugate space nor matches any 
continuum limit. 

6 The covariance of a Gaussian random field that is statistically homoge- 
neous in position space becomes diagonal in the harmonic basis. 
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• The rg_space class comprises all regular Euclidean 
grids of arbitrary dimension and periodic boundary con- 
ditions. Such a grid is described by the number of 
grid points per dimension, the edge lengths of one n- 
dimensional pixel and a few flags specifying the origin of 
ordinates, internal symmetry, and basis type; i.e., whether 
the grid represents a position or Fourier basis. The 
conjugate space of a rg_space is another rg_space 
that is obtained by a fast Fourier transformation of the 
position basis yielding a Fourier basis or vice versa by 
an inverse fast Fourier transformation. 

• The spherical harmonics basis is represented by the 
lm_space class which is defined by the maximum 
of the angular and azimuthal quantum numbers, and 
m, where m max < max and equality is the default. It 
serves as the harmonic basis for the instance of both the 
gl_space and the hp_space class. 

• The gl_space class describes a Gauss -Legendre grid 
on an S 2 sphere, where the pixels are centered at the roots 
of Gauss-Legendre polynomials. A grid representation is 
defined by the number of latitudinal and longitudinal bins, 
ni a t and n\ on- 

• The hierarchical equal area isolatitude pixelization of an 
S 2 sphere (abbreviated as HEALPdQ) is represented by 
the hp_space class. The grid is characterized by twelve 
basis pixels and the n S id e parameter that specifies how 
often each of them is quartered. 

• The nested_space class is designed to comprise all 
possible product spaces constructed out of those de- 
scribed above. Therefore, it is defined by an ordered list 
of space instances that are meant to be multiplied by 
an outer product. Conjugation of this space is conducted 
separately for each subspace. 

For example, a 2D regular grid can be cast to a nesting of 
two ID regular grids that would then allow for separate 
Fourier transformations along one of the two axes. 

B. Fields 

The second fundamental NIFTY class is the field class 
whose purpose is to represent discretized fields. Each field 
instance has not only a property referencing an array of field 
values, but also domain and target properties. The domain 
needs to be stated during initialization to clarify in which space 
the field is defined. Optionally, one can specify a target space 
as codomain for transformations; by default the conjugate 
space of the domain is used as the target space. 

In this way, a field is not only implemented as a simple 
array, but as a class instance carrying an array of values 
and information about the geometry of its domain. Calling 
field methods then invokes the appropriate methods of the 
respective space without any additional input from the user. 
For example, the scalar product, computed by field, dot, 
applies the correct weighting with volume factors as addressed 



fields are defined on different but conjugate domains]^] The 
same is true for all other methods applicable to fields; see 
Tab. [n| for a selection of those instance methods. 

Furthermore, NIFTY overloads standard operations for 
fields in order to support a transparent implementation of 
algorithms. Thus, it is possible to combine field instances 
by -f ,—,*,/,.. . and to apply trigonometric, exponential, and 
logarithmic functions componentwise to fields in their current 
domain. 

C. Operators 

Up to this point, we abstracted fields and their domains 
leaving us with a toolkit capable of performing normalizations, 
field-field operations, and harmonic transformations. Now, we 
introduce the generic operator class from which other, 
concrete operators can be derived. 

In order to have a blueprint for operators capable of han- 
dling fields, any application of operators is split into a general 
and a concrete part. The general part comprises the correct 
involvement of normalizations and transformations, necessary 
for any operator type, while the concrete part is unique for 
each operator subclass. In analogy to the field class, any 
operator instance has a set of properties that specify its domain 
and target as well as some additional flags. 

For example, the application of an operator A to a field s is 
coded as A ( s ) , or equivalently A . times ( s ) . The instance 
method times then invokes _brief ing, _multiply and 
_debrief ing consecutively. The briefing and debriefing are 
generic methods in which in- and output are checked; e.g., 
the input field might be transformed automatically during the 
briefing to match the operators domain. The _multiply 
method, being the concrete part, is the only contribution coded 
by the user. This can be done both explicitly by multiplication 
with a complete matrix or implicitly by a computer routine. 

There are a number of basic operators that often appear 
in inference algorithms and are therefore preimplemented in 
NIFTY. An overview of preimplemented derivatives of the 
operator class can be found in Tab. [TTT1 

D. Operator Probing 

Furthermore, the NIFTY library features a generic 
probing class. The basic idea of probing l23l is to approxi- 
mate quantities of implicit operators that are only accessible at 
a high computational expense by using sample averages. The 
sample pool is generated by the application of random fields 
that are constructed in order to project the quantity of interest. 
For example, an approximation of the trace or diagonal of an 
operator A can be obtained by 



tr[A]«<£U* 



«} 



diag[A] « (£*A£) tt} , 



(11) 
(12) 



where (•){£} I s me sam pl e average of a sample of random 
in Sec. |II-B| and performs basis transformations if the two fields £ with the property (£, P £,q)^y —> S pq for |{£}| — >> oo and 



7 HEALPix homepage 

http : //source forge . net /pro jects/healpix/ 



8 Since the scalar product by discrete summation approximates the integra- 
tion in its continuum limit, it does not matter in which basis it is computed. 



5 



TABLE II 

Selection of instance methods of the NIFTy field class. 



method name 



description 



cast_domain alters the field's domain without altering the field values or the codomain. 

conjugate complex conjugates the field values. 

dot applies the scalar product between two fields, returns a scalar. 

tensor_dot applies a tensor product between two fields, returns a field defined in the product space. 

pseudo_dot applies a scalar product between two fields on a certain subspace of a product space, returns a scalar or a field, depending on the 
subspace. 

dim returns the dimensionality of the field, 

norm returns the L 2 -norm of the field, 

plot draws a figure illustrating the field. 

set_target alters the field's codomain without altering the domain or the field values. 

set_val alters the field values without altering the domain or codomain. 

smooth smoothes the field values in position space by convolution with a Gaussian kernel, 

transform applies a transformation from the field's domain to some codomain. 

weight multiplies the field with the grid's volume factors (to a given power), 
(and more) 



TABLE III 

Overview of derivatives of the NIFTy operator class. 



NIFTy subclass 



operator 

°->- diagonal_operator 
°->- power_operator 

°->- pro jection_operator 

vecvec_operator 
°->- response_operator 



description 



representing diagonal matrices in a specified space. 

representing covariance matrices that are defined by a power spectrum of a statistically homogeneous and isotropic 
random field. 

representing projections onto subsets of the basis of a specified space, 
representing matrices of the form A = aa^ , where a is a field. 

representing an exemplary response that includes a convolution, masking and projection. 



* denotes componentwise multiplication. Since the residual 
error of the approximation decreases with the number of 
samples, one has to find a tradeoff between numerical accuracy 
and computational cost. 

The NIFTY probing class allows for the implementation 
of arbitrary probing schemes. Because each sample can be 
computed independently, all probing operations take advantage 
of parallel processing for reasons of efficiency, by default. 
There are two derivatives of the probing class implemented in 
NIFTY, the trace_probing and diagonal_probing 
subclasses, which enable the probing of traces and diagonals 
of operators, respectively. 

An extension to improve the probing of continuous op- 
erators by exploiting their internal correlation structure as 
suggested in Ref. [ 22 1 is planned for a future version of 
NIFTY. 

IV. Demonstration 

An established and widely used inference algorithm is the 
Wiener filter Q whose implementation in NIFTY shall serve 
as a demonstration example. 

The underlying inference problem is the reconstruction of 
a signal, s, from a data set, d, that is the outcome of a 
measurement process ([2]), where the signal response, R s, is 
linear in the signal and the noise, n, is additive. The statistical 
properties of signal and noise are both assumed to be Gaussian, 



s -r^ Q(s, S) oc exp (— \ 
n ^ Q(n,N). 



5 t<T 



(13) 
(14) 



Here, the signal and noise covariances, S and N, are known 
a priori. The a posteriori solution for this inference problem 



can be found in the expectation value for the signal m = 
( s )( s \d) weight by the posterior P(s\d) . This map can be 
calculated with the Wiener filter equation, 



m 



^y 1 



(15) 



which is linear in the data. In the IFT framework, this scenario 
corresponds to a free theory as discussed in Ref. lfT2lL where 
a derivation of Eq. fl5] ) can be found. In analogy to quantum 
field theory, the posterior covariance, D, is referred to as the 
information propagator and the data dependent term, j, as the 
information source. 

The NIFTY based implementation is given in App. c]^ This 
implementation is not only easily readable, but it also solves 
for m regardless of the chosen signal space; i.e., regardless 
of the underlying grid and its resolution. The functionality of 
the code for different signal spaces is illustrated in Fig. [T] 

The confidence in the quality of the reconstruction can be 
expressed in terms of a la-confidence interval that is related 
to the diagonal of D as follows, 

<r< m > = y/dmg[D]. (16) 

The operator D defined in Eq. fj"5] ) may involve inversions 
in different bases and thus is accessible explicitely only with 
major computational effort. However, its diagonal can be 
approximated efficiently by applying operator probing ( fT2| ). 
Fig. [2] illustrates the ID reconstruction results in order to 
visualize the estimates obtained with probing and to emphasize 
the importance of a posteriori uncertainties. 

9 The Wiener filter demonstration is also part of the NIFTy package; see 
nif ty /demos /demo_excaliwir .py for an extended version. 
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Fig. 1. Illustration of the Wiener filter code example showing (left to right) a Gaussian random signal, the data including noise, and the reconstructed map. 
The additive Gaussian white noise has a variance <j\ that sets a signal-to-noise ratio (c s )^ /a n of roughly 2. The same code has been applied to three 
different spaces, namely (top to bottom) a ID regular grid with 512 pixels, a 2D regular grid with 256 x 256 pixels, and a HEALPix grid with n s id e = 128 
corresponding to 196, 608 pixels on the S 2 sphere. (All figures have been created by NIFTY using the field, plot method.) 



The Wiener filter code example given in App. [C] can easily 
be modified to handle more complex inference problems. 
In Fig. [3] this is demonstrated for the image reconstruction 
problem of the classic "Moon Surface" imag^] During the 
data generation ([2]), the signal is convolved with a Gaussian 
kernel, multiplied with some structured mask, and finally, con- 
taminated by inhomogeneous Gaussian noise. Despite these 
complications, the Wiener filter is able to recover most of the 
original signal field. 

NIFTY can also be applied to non-linear inference prob- 
lems, as has been demonstrated in the reconstruction of 
log-normal fields with a priori unknown covariance and 
spectral smoothness l24l . Further applications reconstructing 
three-dimensional maps from column densities l25l and non- 
Gaussianity parameters from the cosmic microwave back- 
ground l26l are currently in preparation. 



10 Source taken from the USC-SIPI image database at 

http ://sipi.usc. edu/ database/ 



V. Conclusions & Summary 

The NIFTY library enables the programming of grid and 
resolution independent algorithms. In particular for signal 
inference algorithms, where a continuous signal field is to be 
recovered, this freedom is desirable. This is achieved with an 
object-oriented infrastructure that comprises, among others, 
abstract classes for spaces, fields, and operators. NIFTY 
supports a consistent discretization scheme that preserves the 
continuum limit. Proper normalizations are applied automati- 
cally, which makes considerations by the user concerning this 
matter (almost) superfluous. NIFTY offers a straightforward 
transition from formulas to implemented algorithms thereby 
speeding up the development cycle. Inference algorithms that 
have been coded using NIFTY are reusable for similar infer- 
ence problems even though the underlying geometrical space 
may differ. 

The application areas of NIFTY are widespread and include 
inference algorithms derived within both information field 
theory and other frameworks. The successful application of 
a Wiener filter to non-trivial inference problems illustrates the 
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coordinate coordinate coordinate 



Fig. 2. Illustration of the ID reconstruction results. The left panel summarizes the results from Fig. [T] by showing the original signal (red dashed line), the 
reconstructed map (green solid line), and the lcr-confidence interval (gray contour) obtained from the square root of the diagonal of the posterior covariance 
D that has been computed using probing; cf. Eq. fl2) . The middle panel shows the ID data set from Fig.fflwith a blinded region in the interval [0.5, 0.7]. The 
right panel shows again the original signal (red, dashed line), the map reconstructed from the partially blinded data (green solid line), and the corresponding 
lcr-interval (gray contour) which is significantly enlarged in the blinded region indicating the uncertainty of the interpolation therein. 



signal oata reconstructed map 




f 



(centered) convolution kernel mask noise standard deviation 




Fig. 3. Application of a Wiener filter to the classic "Moon Surface" image on a 2D regular grid with 256 x 256 pixels showing (top, left to right) the 
original "Moon Surface" signal, the data including noise, and the reconstructed map. The response operator involves a convolution with a Gaussian kernel 
(bottom, left) and a masking (bottom, middle). The additive noise is Gaussian white noise with an inhomogeneous standard deviation (bottom, right) that 
approximates an overall signal-to-noise ratio (<7 s ) n / (cr n ) Q of roughly 1. (All figures have been created by NIFTY using the field. plot method.) 



flexibility of NIFTY. The very same code runs successfully 
whether the signal domain is an n-dimensional regular or a 
spherical grid. Moreover, NIFTY has already been applied to 
the reconstruction of Gaussian and log-normal fields [24 1. 

The NIFTY source code and online documentation is pub- 
licly available on the project homepage FM 



11 NIFTY homepage 

http : / /www . mpa-garching . mpg . de/ if t /nifty/ 
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Appendix 

A. Remark On Matrices 

The discretization of an operator that is defined on a 
continuum is a necessity for its computational implementation 



and is analogous to the discretization of fields; cf. Sec. II-B 



However, the involvement of volume weights can cause some 
confusion concerning the interpretation of the corresponding 
matrix elements. For example, the discretization of the contin- 
uous identity operator, which equals a ^-distribution S(x — y), 
yields a weighted Kronecker-Delta S pq , 



id = S(x-y)^((5{x-y)) ) = 



V a 



(17) 



where x G Q, p and y G fl q . Say a field £ is drawn from a 
zero-mean Gaussian with a covariance that equals the identity, 
£/(£,id). The intuitive assumption that the field values of £ 
have a variance of 1 is not true. The variance is given by 



(18) 



and scales with the inverse of the volume V q . Moreover, the 
identity operator is the result of the multiplication of any 
operator with its inverse, id = A -1 A. It is trivial to show 
that, if A(x, y) ^ A pq ™ A ^ A ~ x 
A maps as follows, 



an( l J2 q A p qA qr 



S pr , the inverse of 



((A- 1 (x-y)) a ) a ={A- 1 ) pq = 



A- 1 



(19) 



where A~ q x in comparison to (A~ 1 ) pq is inversely weighted 
with the volumes V p and V q . 

Since all those weightings are implemented in NIFTY, users 
need to concern themself with these subtleties only if they 
intend to extend the functionality of NIFTY. 

B. Libraries 

NIFTY depends on a number of other libraries which are 
listed here for completeness and in order to give credit to the 
authors. 

. gfftEO 

. HEALP^ F] (a nd HEALPix (771) 

. LibSHARF|^J[28] (wrapped in PYTHON) 

. NumPy and SciP^] g9) 

• (several other Python standard libraries) 

12 GFFT homepage https://github.com/mrbell/gfft 
13 HEALPy homepage https : / /github . com/healpy /healpy 
14 LiB SHARP homepage 
http : //source forge . net /pro jects/ lib sharp/ 

15 NumPy and SciPy homepage http://numpy.scipy.org/ 
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C. Wiener Filter Code Example 

from nifty import * 

from scipy . sparse . linalg import LinearOperator as lo 
from scipy . sparse . linalg import eg 



class propagator (operator) : 



# define propagator class 



_matvec = (lambda self, x: self . inverse_times (x) . val . flatten () ) 

def _multiply ( self , x) : 

# some numerical invertion technique; here, conjugate gradient 

A = lo ( shape=tuple ( self . dim ()) , matvec=self ._matvec) 

b = x . val . flatten ( ) 

x_, info = eg (A, b, M=None) 

return x_ 

def _inverse_multiply ( self , x) : 
S, N, R = self. para 

return S . inverse_t imes (x) + R . ad j oint_t imes (N . inverse_t imes (R. times (x) ) ) 



# some signal space; e.g., a one-dimensional regular grid 
s_space = rg_space (512, zerocenter=False , dist=0.002) 

# or rg_space ( [25 6, 256]) 

# or hp_space ( 128 ) 



# define signal space 



k_space = s_space . get_codomain ( ) 

kindex, rho = k_space . get_power_index ( irreducible=True ) 



# get conjugate space 



# some power spectrum 

power = [42 / (kk +1) ** 3 for kk in kindex] 



S = power_operator (k_space, spec=power) 
s = S . get_random_f ield (domain=s_space ) 



# define signal covariance 

# generate signal 



R = response_operator ( s_space, sigma=0.0, mask=1.0, assign=None) 
d_space = R. target 



# define response 

# get data space 



# some noise variance; e.g., 1 

N = diagonal_operator (d_space, diag=l, bare=True) 
n = N . get_random_f ield ( domain=d_space ) 



# define noise covariance 

# generate noise 



d = R(s) + n 



# compute data 



j = R . ad j oint_t imes (N . inverse_t imes (d) ) 

D = propagator ( s_space, sym=True, imp=True, para= [ S , N, R] ) 



# define source 

# define propagator 



m = D(j) 

s .plot ( t it le=" signal" ) 
d . cast_domain ( s_space ) 

d . plot (title="data" , vmin=s . val .min ( ) , vmax=s . val . max ( ) ) 



# reconstruct map 

# plot signal 

# plot data 



m . plot (title=" reconstructed map", vmin=s . val . min ( ) , vmax=s . val . max ( ) ) # plot map 



